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We propose a coherent method for the detection and reconstruction of gravitational wave signals 
with a network of interferometric detectors. The method is derived using the likelihood functional 
for unknown signal waveforms. In the standard approach, the global maximum of the likelihood over 
the space of waveforms is used as the detection statistic. We identify a problem with this approach. 
In the case of an aligned pair of detectors, the detection statistic depends on the cross-correlation 
between the detectors as expected, but this dependence dissappears even for infinitesimally small 
misalignments. We solve the problem by applying constraints on the likelihood functional and obtain 
a new class of statistics. The resulting method can be applied to data from a network consisting of 
any number of detectors with arbitrary detector orientations. The method allows us reconstruction 
of the source coordinates and the waveforms of two polarization components of a gravitational wave. 
We study the performance of the method with numerical simulation and find the reconstruction of 
the source coordinates to be more accurate than in the standard approach. 
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I. INTRODUCTION 

Several gravitational wave (GW) detectors are now op- 
erating around the world, including both laser interfer- 
ometers PISHE! an( A resonant mass detectors 0- Com- 
bining the data from such a network of detectors can ben- 
efit both the detection of GW signals and estimation of 
signal parameters. Unlike real GW signals that would oc- 
cur in coincidence across all detectors in a network, most 
background events due to instrumental and terrestrial 
disturbances are expected to be local to each detector 
and, therefore, can be rejected by analysing data from a 
network. Given a network with different orientations and 
locations of the detectors, GW sources can be localized 
on the sky and the waveforms of the two independent 
GW polarization components can be reconstructed. 

Methods for the analysis of data from a network of GW 
detectors can be divided into two classes: coincidence and 
coherent methods. In coincidence methods, first, a search 
for GW signals is carried out for individual detectors and 
a list of candidate events is generated. Then a subset of 
events is selected by requiring temporal coincidence of 
events between the detectors. In coherent methods, one, 
first, combines the detector responses and then analyzes 
the combined data to generate a single list of events. 

Networks of detectors are particularly important for 
searches of gravitational wave burst signals. These are 
defined to be broadband signals that may come either 
from unanticipated sources or from sources for which no 
reliable theoretical prediction exists for signal waveforms. 
Potential astrophysical sources of burst signals are stellar 
core collapse in Supernovae |gj, mergers of binary neutron 
star or black hole systems Q and Gamma Ray Burst 
progenitors [8J- 

The first coherent method for burst searches with a 
network of three misaligned detectors was proposed by 
Gursel and Tinto |9(. In this method, the detector re- 



sponses are combined into a functional, which attains its 
minimum at the correct direction to the source. The min- 
imization of the functional allows one to reconstruct the 
source coordinates and two polarization waveforms of the 
burst signal. 

Flannagan and Hughes IToll considered maximization 
of the likelihood functional [111 [12J as a means of recon- 
structing source direction and polarization waveforms. 
Anderson et al [Rj extended this approach to derive a 
detection statistic called excess power. It is obtained 
by integrating the likelihood functional, weighted by a 
Bayesian prior probability density, over the space of all 
waveforms. In this paper, we refer to signal detection 
and reconstruction based on the global maximum of the 
unweighted likelihood functional as the standard likeli- 
hood method. Another coherent method, proposed by 
Sylvestre [3, starts with the ad hoc approach of form- 
ing a linear combination of data from a network of detec- 
tors. The combination coefficients are then adjusted to 
construct a quadratic detection algorithm that satisfies 
certain well defined criteria. 

Arnaud et al |15| have numerically explored the issue of 
statistical performance of coherent and coincidence meth- 
ods and find that the former are more efficient than the 
latter for burst si gnal s. In the case of signals with known 
waveforms, Finn |l6| has shown using simulations that 
a coherent method can also be robust when confronted 
with non-Gaussian noise. 

Coherent methods can also be used for rejecting coin- 
cident background signals. Cadonati has proposed 
a cross-correlation test, called r- statistic, for pairs of 
aligned detectors as a follow up consistency check on a 
coincidence analysis ^[ . Rakhmanov and Klimenko frof 
have proposed the mixed correlations method that ex- 
tends the cross-correlation test to a network of three or 
more misaligned detectors. Wen and Schutz [2(j have re- 
cently generalized the coherent approach of Gursel and 
Tinto [9fl to create a method for rejecting background 
coincident signals with a network of arbitrary detectors. 

In this paper, we propose a method for the coherent de- 
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tection and reconstruction of burst si gnal s that is based 
on the use of the likelihood ratio flU Il2j. Our analysis 
differs from m an important way. We identify and 

solve a problem with the standard likelihood analysis, 
first spotted in [2l|. The problem, which we call the two 
detector paradox, is that the maximum likelihood ratio 
statistic for misaligned detectors does not reduce, con- 
trary to physical intuition, to the statistic for co-aligned 
detectors in the limit of small misalignment angles. The 
latter statistic depends on the cross-correlation of detec- 
tor outputs whereas the former does not. We show that 
the problem originates in the maximization of the likeli- 
hood ratio functional over all signal waveforms including 
those to which a detector network may not actually be 
sensitive. We propose a solution to this problem that is 
based on constraints imposed on the GW signal wave- 
forms. 

The constrained maximization of the likelihood func- 
tional yields new detection and reconstruction methods 
which we call the constraint likelihood methods. Unlike 
the Giirsel and Tinto method, the constraint likelihood 
methods can be used for arbitrary networks, including 
networks consisting of two detectors. The performance 
of these methods is studied in comparison with the stan- 
dard likelihood method by using the numerical simula- 
tions. In the simulation we use networks of interferomet- 
ric detectors consisting of LIGO 4 km detector in Hanford 
(HI), LIGO 4 km detector in Livingston (LI), GEO-600 
detector (Gl), TAMA detector (Tl) and VIRGO detec- 
tor (VI). We find that the constraints employed in this 
paper enhance the detection efficiency of the likelihood 
method. For detected sources, the constraints signifi- 
cantly improve accuracy of the source localization. 

The rest of the paper is organized as follows. ScctionlTTI 
lays out much of the basic notation and conventions used 
in the paper. In Section lllll we provide an overview of 
the standard likelihood approach and its application to 
burst signals. Section|^describes the two detector para- 
dox that appears in the standard likelihood approach. 
The origin of this problem is discussed in Section In 
Section WH we derive the constraint likelihood methods. 
The results from numerical studies of the performance of 
these methods are described in Section IVIII 



II. DETECTOR RESPONSE TO 
GRAVITATIONAL WAVES 

A. Gravitational wave signal 

Gravitational waves are described by a symmetric ten- 
sor of second rank hij(t), which is usually defined in 
the transverse-traceless gauge [S^- It takes particularly 
simple form in the coordinate frame associated with the 
wave. In this coordinate frame (the wave frame), a grav- 
itational wave propagates in the direction of z axis and 
it can be described with the waveforms h+(t) and h x (t) 
representing two independent polarizations components 



of the wave. 

In addition to the waveforms h+(t) and h x {t) we will 
use complex waveforms defined as 

u(t) = h+(t)+ih x (t), (I) 
u(t) = h+(t) - ih x (t). (2) 

In what follows tilde will always denote complex conju- 
gation. The GW waveforms u(t) and u(t) are eigenstates 
of the rotations around z-axis in the wave frame. We 
denote this particular rotation by R z (ip), where ip is the 
rotation angle. The rotation R z (ip) generates equivalent 
waveforms which arc different representations of the same 
gravitational wave. 

We define the sum-square energy [^j carried by the 
gravitational wave as 

/OC />OG 
(h\{t) + h 2 x {t)) dt = / u(t)u(t)dt. (3) 
-OC J — oo 

Note that the sum-square energy is invariant under the 
rotation R z . 

B. Detector response 

The response of the interferometer to an arbitrary 
gravitational wave hij (t) is given by 

m = lT ij h ij (t), (4) 

where is the detector tensor 23] . In the wave frame 
the detector response is a linear superposition of two GW 
polarizations 

C(t) =F + h + (t) + F x h x (t). (5) 

where the coefficients F + and F x are known as antenna 
patterns. 

To calculate the antenna pattern we introduce the 
Earth-centered frame described in [24| . In this frame the 
detector location is defined by a radius-vector r point- 
ing to the detector and its orientation is described by 
two unit vectors a and b along the detector arms. The 
vectors a and b define the detector tensor 

Tlj = a, a,j - bi bj, i, j = 1, 2, 3, (6) 

where the indices correspond to spatial coordinates x, 
y and z respectively. The direction to the GW source 
is defined in the Earth-centered frame by two spherical 
angles (j> (longitude) and 9 (lattitude). The rotational 
transformation which connects the Earth-centered frame 
with the wave frame is given by 

TL(<f>,6) = -R v (6)R z (<p). (7) 

It defines the detector tensor in the wave frame 

T(0,0)=R(0,0)T'R(0,0) T (8) 



3 



Omitting the explicit dependence on the angles, the an- 
tenna patterns corresponding to the h + (t) and h x (t) po- 
larizations are calculated as follows: 



- \{Ti2 + T 21 ). 



(9) 
(10) 



The detector response can be conveniently expressed 
in terms of the complex waveform u: 

£ = Au + Au, (11) 
where A and A are the complex antenna patterns: 

(12) 
(13) 



A = \ (F + + iF x ), 
1 



A = - (F+-iF x ). 



A rotation i? z (V>) in the wave frame induces the trans- 
formation of the detector antenna patterns and the GW 
waveforms 



A' = e 2 ^ A, 
v! = e 2 ^ «, 

but the detector response is invariant under the rotation 



(14) 
(15) 



III. LIKELIHOOD ANALYSIS OF 
GRAVITATIONAL WAVE DATA 

In this section, we present a brief overview of the stan- 
dard likelihood approach to the detection and reconstruc- 
tion of gravitational wave burst signals using a network 
of detectors. Though the scientific content of this sec- 
tion is essentially the same as the results in 0, , our 
derivation and the notation we use are rather different. 
These will aid in a clearer exposition of our main results 
in subsequent sections. The reader is referred to for 
a textbook level discussion of the statistical theory of 
signal detection used in this paper. 



A. Overview 

Consider an observable that is a finite data segment 
x = {x[l], x[2], . . . , x[AT]} from a noisy time series. The 
simplest detection problem is to define a decision rule for 
selecting one of two mutually exclusive hypotheses, Hq 
(null hypothesis) or Hi (alternative hypothesis), about 
the data x. Under the Hq and Hi, x is a realization 
of a stohastic process described by the joint probability 
density p(x | Hq) and p(n\Hi) respectively. 

Any decision rule will incur two types of errors: false 
alarm - Hi is selected when Hq is true, and false dismissal 
- Hq is selected when Hi is true. Each error will have 
a probability associated with it, namely the false alarm 



and the false dismissal probabilities Qq and Qi respec- 
tively. In order to select the best decision rule, several 
criteria have been proposed out of which the Neyman- 
Pearson criterion is the most suitable for detection of 
gravitational waves. According to this criterion, the op- 
timal decision rule has the least Qi f° r fixed Qq. The rule 
accepts Hi (Hq) when the likelihood ratio, A(x), defined 
as 



A(x) = 



P(x|#i) 
p(x|ifo) 



(16) 



is greater (less) than a threshold value that is fixed by 
the specified Q . 

In the case of the GW data analysis, Hq is the hypoth- 
esis "a GW signal is absent" and Hi is "the GW signal £ 
is present" . For a stationary, Gaussian white noise with 
zero mean the corresponding joint probability densities 
are 



N 1 / 2ri 



i=l 
N 



27TCT 
1 



(17) 



i—l x / 

where a is the standard deviation of the noise. The log- 
arithm of the likelihood ratio can be expressed as 



N 1 / 1 
£ = ln(A(x)) = ^- i «i]--e 2 W 



(19) 



In the rest of the paper, we will be concerned only with 
C which will be referred to as simply the likelihood. 

The situation with two mutually exclusive hypotheses, 
outlined above, is the simplest one. In general, as in the 
case of GW analysis, the observed data x can be a re- 
alization of one among several joint probability densities 
p(x\Hi), i = 0,1,2,..., where, as usual, Hq is the null 
hypotheses and Hi are the alternative hypotheses. Cor- 
respondingly, the probabilities for false alarm and false 
dismissal can be assigned but now the false dismissal 
probabilities Qi are hypothesis specific. 

One possible generalization of the Neyman-Pearson 
criterion could be to select that decision rule which min- 
imizes all probabilities Qi for a fixed false alarm prob- 
ability Qq. It turns out that, in general, no such rule 
is possible llj. Another approach is to generalize the 
likelihood ratio test itself by constructing a functional 



A, 



p{x\Hi) 



p(x|#o) 



(20) 



and comparing it with a threshold. This test, called the 
maximum likelihood ratio (MLR) test, tends to outper- 
form any other ad hoc test. However, it is important to 
note that the MLR test itself does not have a formal proof 
of optimality. Therefore, it is possible that modifications 
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of the MLR test, as presented in this paper, can lead to 
better performance. 

One of the applications of the MLR test is the detec- 
tion of gravitational waves from the inspiral of compact 
binaries |16l |25| . In principle, the waveforms of the GW 
signals can be calculated to arbitrary precision given the 
parameters of the binary system. The set of alternative 
hypotheses now becomes a continuum that is identified 
with the space of binary parameters. The likelihood ratio 
A(x|ilj) can, therefore, be expressed as a function over 
the binary parameters. The MLR statistic is obtained 
by maximizing the likelihood ratio over these parame- 
ters and reaches its maximum for the best match of the 
corresponding waveform to the data. 

In contrast to binary inspiral signals, where the num- 
ber of parameters is small, the parameters characterizing 
burst signals are essentially the signal amplitudes them- 
selves at each instant of time. Thus, for burst signals, 
the number of parameters can be very large. Formally, 
however, the concept of the likelihood ratio can still be 
used for burst signals. In this case, the likelihood ratio 
is A(x|£), where £ is the detector response to the burst 
signal. The application of the MLR test to burst sig- 
nals involves maximization over each sample inde- 
pendently 



We also define the network output time series X which 
combines the output time series X& from individual de- 
tectors 



fe=i 



err 



(25) 



To simplify equations, we will replace summation over 
any sampled time series s[i] with (s). With these new 
notations the likelihood functional can be written as 



C = < uX + uX — g r uu 



g c u 2 + g c u 2 



(26) 



where the X and g c are complex conjugates of X and g c 
respectively. 



C. Solution for GW waveforms 

The equations for the GW waveforms are obtained by 
variation of the likelihood functional: 



sc 

Su 



SC 

6 u 



0, 



(27) 



B. Network likelihood 



which results in two linear equations for u and u 



So far, we have considered a time series x at the output 
of a single GW detector. The entire formalism outlined 
above can be extended to a network of detectors. Let the 
data from the k th detector be x^ = {xfc[l], £fe[2], . . .} and 
the detector response to the gravitational wave be 



£fc[£] = u[i]A k + u[i]A k 



(21) 



We will assume that the noise in different detectors is 
independent. Then the likelihood ratio becomes, 



K N 



=1 u k 



(22) 



where K is the number of detectors in the network. For 
detectors illuminated by the same GW source, the detec- 
tor responses are not independent. Therefore, the vari- 
ation of the likelihood functional is performed over the 
sampled amplitudes u[i] and u[i\. 

To characterize the angular and strain sensitivity of 
the network, we introduce the network antenna patterns 



K A A K 

k=l k k=l 



At 



(23) 



where g r is real and g c is complex. Similarly to the an- 
tenna patterns A k for a single detector, they describe the 
network response to the gravitational wave: 



R(u) = g r u + g c u . 



(24) 



The solution is 



X = g r u + g c u, 
X = g r u + g c u. 



g r X - g c X 



fir 



\f]e 



(28) 
(29) 



(30) 



Note, the solution u s satisfies the condition X = R(u s ), 
where R(u s ) is the network response to the gravitational 
wave (see Ea. (|24|) L Equations 1281 and 1291 can also be 
written in the matrix form 



Re(X) 
lm(X) 



where the matrix Mr is given by 

Mr 



g r + Re(p c ) lm(0 c ) 
Im (.g c ) 9r ~ Re(g c ) 



(31) 



(32) 



D. Maximum likelihood ratio statistic 

The maximum likelihood ratio statistic is obtained by 
substitution of the solution u s in Eq. I|2()l) 



2g r X r — g c X c — g c X c 

H9 2 r-\9c\ 2 ) 



(33) 
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where the quantities X c and X r are defined by 

K K 

X c = ^^ijAj, (34) 

«=1 3=1 
if if 

i=l j=l 

and A c is complex conjugate of X c . The data matrix 
is calculated for the detector output x; and xj scaled by 
the variances of the detector noise 



(Xj (t)Xj(t + Tjj)) 



(36) 



The data matrix depends on the gravitational wave time 
delays between the detectors. The time delays, in 
turn, depend on the coordinates of the source on the 
sky 9 and (p. The diagonal elements of the data matrix 
represent the power terms and the non-diagonal elements 
represent the cross-correlation terms. 

There is a simple geometrical interpretation of the 
MLR statistics. At any instance of time, the GW wave- 
form u and the network output X can be viewed as vec- 
tors u and X in the complex plane. Then the MLR 
statistics, given by Ea. (|33[l . is the inner product 



u s X + u s X) = (u s -X 



(37) 



which is a projection of the solution u s onto the data X. 
Note, that the projection is the estimator of the total 
signal-to-noise ratio of the GW signal detected in the 
network 



K 



SNR tot = J2 



2 (u s • X) 



(38) 



fc=i 



IV. TWO DETECTOR PARADOX 

So far we have described the standard likelihood ap- 
proach for the detection and reconstruction of burst GW 
signals wherein the likelihood ratio is maximized inde- 
pendently over each signal sample. Though attractive 
because both a detection and estimation method are ob- 
tained simultaneously, there is a problem with this ap- 
proach when applied to a network of two detectors. The 
problem, which we call the two-detector paradox, is de- 
scribed in this section. 

Let us consider a network of two detectors in two con- 
figurations: (A) aligned detectors and (M) misaligned 
detectors. The detectors in the configuration A have 
the same antenna patterns. In this case the detector 
responses are the same in both detectors and we con- 
sider the GW signal as the scalar wave £. The likelihood 
functional is then 



L A = 



too , fa® (e) ( i 



o\ 



+ 



o\ 



1 



(39) 



where xi, X2 are the detector outputs and o~\, o~i are the 
standard deviations of the detector noise. The solution 
of the likelihood variation problem is 



1 



1 



(40) 



The MLR statistics for two aligned detectors is obtained 
from Ea. (|39|l by substituting £ with the solution 



L 



(41) 

As expected, the MLR statistic for two aligned detectors 
includes both the power and the cross-correlation terms. 
For two arbitrary misaligned detectors the MLR statistic, 
given by Ea. (|33|l . reduces to 




'<*?> , <*!>' 



(42) 



which includes the power terms only. 

The two detector paradox is that the statistic Lm does 
not include cross-correlation between the detectors even 
for a small misalignment. This is highly counterintuitive 
since one expects that the response of detectors to the 
same GW source will differ only infinitesimally when the 
detectors are infinitesimally misaligned. Hence, as in the 
case of La, one would expect that the cross-correlation 
term will benefit detection and that its importance will 
decline only gradually as the detectors are misaligned. In 
other words, the functional Lm is expected to approach 
La in the limit of perfect alignment. 

The origin of the two detector paradox is easily seen. 
For the aligned case, the standard likelihood ratio ap- 
proach has the prior information that both detector re- 
sponses are identical. Hence, the cross-correlation term 
is guaranteed to have positive mean and, thus, should 
improve the detectability of GW signals. While, for the 
misaligned case, it is always possible to specify two ar- 
bitrary responses and invert them to obtain some h+(t) 
and h x (t) components of the GW signal. The standard 
MLR statistic, therefore, does not benefit from having 
the cross-correlation term since now it can contribute 
pure noise to the statistic. Hence, this term disappears 
from the MLR statistic. The fact that the standard like- 
lihood approach does not exhibit the expected continuity 
for the case of two detectors indicates that this approach 
may not be the best one for a general network of GW 
detectors also. 



V. NETWORK RESPONSE 

To resolve the two detector paradox we take a closer 
look at how the GW signal and the detector noise con- 
tribute to the MLR statistic. In this section we show 
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that the detection of two GW components can be consid- 
ered as two independent measurements equally affected 
by the detector noise but conducted with different an- 
gular and strain sensitivities of the detectors. Being an 
ad hoc method, the maximum likelihood may not be an 
optimal approach in this situation. For example, if the 
network is sensitive only to one signal component (as in 
the case of co-aligned detectors) the measurement of the 
second component does not benefit the GW detection, 
but rather adds noise to the measurement. In the next 
section we propose a solution to the problem and derive 
the detection statistics, which continuously bridge the 
cases of aligned and misaligned detectors. 

A. Network response to gravitational waves 

As we mentioned in Section lllBI the detector response 
is invariant under rotations R z in the wave frame. Con- 
sequently, all measurable quantities, including the likeli- 
hood functional, are invariant as well. We have a freedom 
to select an arbitrary wave frame by applying the rotation 
R z (ip), where ip is the rotation angle. The rotation in- 
duces the transformation of the GW waveforms and the 
detector antenna patterns (see Eq.(JI3J), as well as the 
transformation of the network parameters: X — » Xe %2 ^ 
and g c — > g c e li ^ . In general, the rotation angle ip can 
be selected individually for each instance of time and 
for each point in the sky. By applying the rotation 
R z (— 7/4), where 7 is the phase of g c , we selected a wave 
frame in which both network antenna patterns are real 
and positively defined. We call this particular coordinate 
frame the dominant polarization frame . 

As follows from Eg. (|24fl . for a GW signal u defined in 
the dominant polarization frame, the network response 
is 

R = (9r + \9c\) hi +i(g r - \g e \) hi , (43) 

where hi and hi are the real and imaginary compo- 
nents of the signal. We will distinguish them from the 
GW polarizations h + and h x defined for an arbitrary 
wave frame. Note, the coefficients in front of hi and hi 
are the eigenvalues of the network response matrix Mr 
fEa.l|32[l). which takes a diagonal form in the dominant 
polarization frame 

Mn = g(\ °) • (44) 

The coefficient 

.9 = 9r + \9c\ (45) 

characterizes the network sensitivity to the hi wave. The 
sensitivity to the second component hi is eg, where e is 
the network alignment factor: 



The alignment factor e shows the relative sensitivity of 
the network to the GW components hi and hi. Note 
that < e < 1. The total signal-to-noise ratio of the 
GW signal detected in the network is 

SNR to t =2g((hl)+e(h 2 2 }) , (47) 

where (h\) and (hi) are the sum-square energies carried 
by each component (see Eq.©). Therefore, to be de- 
tected with the same signal-to-noise ratio, the hi wave 
should carry 1/e times more energy than the hi wave. 

Both the network sensitivity and the alignment factor 
depend on the angular and the strain sensitivities of the 
detectors. The alignment factor reflects also the angu- 
lar alignment of the detectors. For co-aligned detectors 
e = and the hi component of the GW signal can not 
be detected. Even for detectors with large angular mis- 
alignment, depending on the sky coordinates 8 and cj>, 
the alignment factor may take small values indicating 
that the detectors are effectively aligned. For example, 
Figure shows the alignment factors as a function of the 
sky coordinates calculated for several network configura- 
tions consisting of the HI, LI, Gl, VI and Tl detectors. 
For simplicity, we assume that the detectors have the 
same strain sensitivity. The example shows, that for the 
closely aligned Hl-Ll detectors, the alignment factor is 
close to zero everywhere, except for a few small patches 
on the sky. The more detectors are added to the net- 
work, the larger is the area on the sky with large values 
of e. But even for the network of five detectors (Hl-Ll- 
G1-V1-T1), the factor e remains small for a considerable 
fraction of the sky area, where the network is much less 
sensitive to the hi wave, than to the hi wave. Assuming 
that both components carry on average the same energy, 
the hi wave is suppressed by the factor of e. Therefore, 
the hi component adds little to the total signal-to-noise 
ratio SNRtot for GW signals originating from areas on 
the sky with small values of e. 

The coefficient g defines the overall sensitivity of the 
network to the gravitational waves. Figure [3 shows the 
network sensitivity calculated as a function of the sky 
coordinates for several network configurations. As we 
expect, adding more detectors reduces the sky area where 
the network is blind to gravitational waves. 

B. Two components of the likelihood functional 

In the dominant polarization frame the likelihood func- 
tional can be written as 

£(u) = ^uA 7 + uA 7 - g r uu ~ ^-(u 2 + u 2 )^ . (48) 

where X 1 = Xe^ %1 ^ 2 . Expressed in terms of hi and hi, 
it can be written as C(hi, hi) = Ci(hi) + Ci(hi): 

Ct = 2 (jX\ cos(/3)/ii - , (49) 

£ 2 = 2(|X|sm(/3)/ l2 -|ftl) , (50) 
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FIG. 1: Alignment factors for the detector networks listed in 
the order from top to bottom: Hl-Ll (upper plot), H1-L1-G1, 
H1-L1-G1-V1, H1-L1-G1-V1-T1 (bottom plot). 



FIG. 2: Sensitivity of the detector networks listed in the order 
from top to bottom: Hl-Ll (upper plot), H1-L1-G1, Hl-Ll- 
Gl-Vl, H1-L1-G1-V1-T1 (bottom plot). 



where \X\ is the amplitude and (3 is the phase of the data 
vector Xn/. The solutions for the h\ and hi are obtained 
by the variation of the C\ and Ci functionals: 

fn = -\X\ cos(/3) , h 2 = —\X\ sin(/3) , (51) 
9 £5 

The MLR statistic can be calculated separately for each 
component 



L\ = -(\Xfcos\p)) = 2Xr + e f Xc ,(52) 

9 4 3 



L^ld^sin^/?))^^- 6 ^-^ 



4eg 



.(53) 



The statistics L\ and Li are the estimators of the signal- 
to-noise ratio of two GW components detected in the 
network. 



C. Detector noise 

The detector output is a sum of the detector noise 
rtfc and the detector response If no GW signal is 
present than the network output is 



K 



x n = j2 



n k A k 



(54) 



k=l 



which follows from the definition of the data vector X 
(see Eq.JSJ)). In this follows from Ea. lj30|l . the 

likelihood variation procedure produces the non-zero so- 
lutions u s for the GW waveforms and the MLR statistic 
is the biased estimator of SNR to t- The reconstructed 
sum-square energy 



En = <ft? n > + (hl n ) = J (L lB + ^) 



(55) 
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is biased as well, where L\ n and L\ n are the MLR statis- 
tics due to the detector noise. The ensemble average E n 
can be easily calculated when the detector noise is white 
and Gaussian. Indeed, in this case the mean of the data 
matrix is 

A~oc^, (56) 

(TiCTj 

where Sij is the Kronecker delta. The average MLR 
statistics due to the detector noise is than 

I^ = I^ocl/2. (57) 

As one can see, in average the detector noise introduces 
the same bias for each signal component. From Eq.J53J) it 
follows that the reconstructed energy in the second com- 
ponent is proportional to 1/e and it diverges when e — > 0. 
Therefore, for small values of e, the likelihood variation 
procedure may result in the un-physical solutions for the 
signal component hi . This is the root of the two detector 
paradox described in Section ITVl 

The statistics L\ and Li can be considered as two in- 
dependent measurements of the GW components hi and 
hi. Indeed, the measurements are uncorrelated, and their 
fluctuations are characterized by the same variances of 
the noise, which follows from the equations 

\L\n — Lin) (Lin — Lin) = 0, L\ n — L\ n . (58) 

When the value of the alignment factor is small, the stan- 
dard MLR statistic L max is not the optimal estimator of 
SNRtot; because the second component adds pure noise 
into the measurement. 



VI. CONSTRAINT LIKELIHOOD 

We have seen that for the standard likelihood approach 
the problem arises when there is a large asymmetry 
(e -C 1) in the detection of two GW components. In this 
case we could find better estimators for the GW wave- 
forms and for the total signal-to-noise ratio of the GW 
signal detected in the network. The construction of such 
estimators depend on our assumptions about the GW 
signals. Mathematically these assumptions can be im- 
plemented as constraints applied to the likelihood func- 
tional. The purpose of the constraints is to exclude the 
un-physical solutions arising from the variation of the 
likelihood functional. By removing such solutions from 
the waveform parameter space we expect to sacrifice a 
small fraction of the real GW signals, while considerably 
improve the detection for the rest of the sources. Below 
we consider examples of the likelihood constraints that 
can be used in the analysis. 

A. Hard constraint 

Given a source population, we could expect that in av- 
erage both signal components hi and hi carry about the 



same energy. For example, for binary sources, the grav- 
itational waves are emitted with the random inclination 
angles. For waves in the dominant polarization frame, 
which is oriented randomly with respect to the source 
frame, the ensemble mean of the sum-square energies of 
two components satisfies 



(hi) = (hi). (59) 

For areas in the sky where the network alignment fac- 
tor is small, for most of the sources the detected energy 
will be dominated by the first component (see Eg . (|47() ) . 
For example, for a network consisting of three interfer- 
ometric detectors HI, LI and Gl the alignment factor 
is less then 0.1 for approximately 50% of the sky area. 
Therefore, the noisy component hi can be entirely ig- 
nored for those sky locations where e is less then some 
threshold eo- This requirement impose a constraint on 
the reconstructed GW waveforms and, therefore, on the 
MLR statistic. For a given sky location we define the 
hard MLR statistic as 

Lw = {/ X It?' ( 6 °) 

When the threshold e = 1, the MLR statistics is defined 
by the first signal component only. In the limit of a 
small alignment angle between the detectors (e — > 0), 
the hard constraint statistics converges to the statistics 
for co-aligned detectors thus resolving the two detector 
paradox. 

The hard constraint is a good approximation in the 
case of closely aligned detectors, such as the network of 
the Hl-Ll detectors. The simulation results (see Sec- 
tion [^nj show, that the L\ is a reasonably good statis- 
tic even for a network of the H1-L1-G1 detectors with 
large angular misalignment between the LIGO and GEO 
detectors. However, if the detection statistic Li is used, 
the search algorithm is entirely inefficient to a GW signal 
when h\ = 0. Although, such GW signals are quite un- 
likely (due to random relative orientations of the source 
and the dominant polarization frames) , and for small val- 
ues of e they may not be detected anyway (unless the hi 
component is very strong), in the next section we intro- 
duce a different constraint, which is free from this prob- 
lem. 



B. Soft constraint 

As we mentioned in Section IV Bl the unconstrained 
MLR statistic L max is a sum of the statistics Li and L2, 
which can be written as 

L max = ^(\X\ 2 (l + 5)), <5=l-isin 2 (/3). (61) 

If the detector output is dominated by the GW signal, 
and assuming that both GW components carry about the 
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same energy, the ensemble average for sin 2 (/3) is 



sin 2 (/3)« £ 7(l + ^), 



(62) 



which follows from the expression for the network re- 
sponse (see Ea.l|4Hjl). It means that in average 



<5« e- 



1 - e 

IT? 



(63) 



and the second term in Ea. (|61|l is much less then 1. On 
contrary, for the detector noise 



sin' 



and respectively 



'(/?) 



1 - e. 



(64) 



(65) 



Therefore, the noisy second term in Ea. l|61(l can be omit- 
ted, resulting in the statistic, which we call the soft MLR 
statistic 



L soft = -{\X\ 2 ) = L 1 +eL 2 
9 



(66) 



There is a simple statistical justification of this result. 
Since the statistics L\ and L 2 are two uncorrelated Gaus- 
sian random variables with the mean \x\ and 1x2, and 
the variance v, the joint probability P(L±, L2, /Zi, fi2, v) 
belongs to the Rayleigh distribution. For the assump- 
tion above (see Eq.lJSSJ), we expect that \i\ = gE and 
/i2 = egE, where E is the GW sum-square energy. Then 
the best estimator for SNR to t is obtained by maximizing 
P over E, which gives the statistic L so f t . 

To obtain the solution for the GW waveforms, one 
should impose a constraint on the likelihood functional 
itself. The constraint can be integrated into the varia- 
tion procedure by the method of the Lagrange multiplier 
|27| . In this method, first, we have to obtain the con- 
straint equation. The soft constraint can be constructed 
by requiring that 



9 (hl)+eg(hl)=0 ) -i(\X\ 2 ) 



(67) 



which limits the sum-square energies (h^) and 
Since, the constraint is applied to the /12 component only, 
we can replace the hi with the solution for the first com- 
ponent and re-write the constraint as 



e .9<^)--<m 2 sin 2 (/3))=0 

g 



(68) 



The solution for the second GW component h,2 is triv- 
ially obtained by the constraint variation of the likelihood 
functional £ 2 



1 



tasoft = sin(/3) 
VC5 



(69) 



As one can see, the constrained solution is the standard 
solution h>2 , multiplied by a penalty factor of ^fe, which 
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FIG. 3: Sky maps of the likelihood statistics L ma x (top), Lhard 
(middle) and L so ft (bottom) for the detector network Hl-Ll- 
Gl. The injected signal SNR is 17 (HI), 20 (LI) and 9 (Gl). 
The source is located at 6 = 50° and d> = 280° . 



reduces both the noise and the signal contribution from 
the second component to the MLR statistic at small val- 
ues of e. Obviously, it reduces the sensitivity of the L so f t 
statistic to a particular class of GW signals with hi — 0, 
described in Section IVI Al For these signals, to be de- 
tected with the same false alarm rate, the £ so ft statistic 
requires (l + e 2 )/2e times more powerful GW signal, then 
the standard L max statistic. For example, for e = 0.1 
the degradation of the strain sensitivity is by a factor of 
2. But it happens only for a small fraction of the GW 
sources. Compare to the standard likelihood method, for 
most of the sources we expect to improve the detection 
sensitivity if the I/ so ft statistic is used. 



C. Network sky maps 

In un-triggered burst searches the coordinates of the 
source, 9 and <f> are free parameters. In this case, the 
detector responses, the likelihood statistics and the re- 
constructed waveforms become functions of 9 and 6 or 
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skymaps. For example, the skymaps corresponding to 
different statistics are shown in Figure|31 (For details see 
Section ED). 

For a given location in the sky, the value of the likeli- 
hood statistic indicates how consistent the data is with 
the hypothesis that a GW signal originates from that lo- 
cation. The coordinates 6 and <f> which yield maximum 
for the likelihood statistic correspond to the most prob- 
able location of the source. The maximum value of the 
statistic is then used for detection. By setting a threshold 
on the maximum likelihood value one can decide on the 
presence or absence of a gravitational wave signal in the 
data as described in Section ITTT1 Given the most prob- 
able source coordinates the waveforms h+(t) and h x (t) 
are reconstructed as described in Section IvTl 



VII. NUMERICAL SIMULATION 

We have outlined a general method for using the 
MLR statistic in conjunction with constraints limiting 
the space of the GW waveforms. The method is intended 
for application to burst searches with networks of gravi- 
tational wave detectors. The performance of the method 
and the effect of the constraints can be analyzed us- 
ing numerical simulations with modeled waveforms. The 
present simulation is similar to the one previously used 
for estimating the performance of the mixed correlation 
method [l9j . 



A. Simulation procedure 

The two-polarization waveforms which represent burst 
gravitational waves used in the simulation are taken from 
the numerical models of the merger phase of coalescing 
binary black holes (BH) (2^. These waveforms form a 
one-parameter family BH-A/, where M is the total mass 
of the binary system in units of solar mass. The results 
below correspond to M = 100. In the simulations we 
generate the detector noise which is Gaussian and white. 
The variance of the noise is selected to be the same for 
all detectors. 

A typical simulated data segment has the duration of 1 
second and consists of N — 4096 data samples. For calcu- 
lation of the data matrix (see Ea. i|rSfi|l 'l we set the integra- 
tion window of 85 ms, which is substantially greater than 
the duration of the signal. The magnitude of the simu- 
lated signals is controlled by the overall gain G, which is 
varied from to 10, whereas the magnitude of the noise 
is kept fixed. 

Due to different orientation of the detectors with re- 
spect to the incoming gravitational wave, the detectors 
responses are different (see Eg. Hill) ). To characterize the 
magnitude of the signal in any given detector we define 



the signal-to-noise ratio: 

SNR = r W$L d f _> J_ ^ ^ (ti ), (to) 

where S(f) is the power spectral density of the noise. 
For white Gaussian noise S(f) = a 2 / f s , where f s is the 
sampling rate. 

For any given detector in the network, the magnitude 
of the signal varies significantly depending on the source 
location in the sky. We therefore choose the location 
of the simulated sources at random, with a uniform dis- 
tribution over the sky. We also choose the polarization 
angle ip at random from the interval [0°,360°]. With 
these choices the simulation gives us an estimate of the 
performance of the detection algorithms without the bias 
which can be introduced by the particular choice for the 
source location or its polarization angle. 

The simulation consists of series of tests corresponding 
to different values of SNR (controlled by G). For each 
value of G a total number of 10,000 injections were made. 
To characterize the strength of the signal in each detector 
for the entire test we introduce the sky-average SNR, 
denoted by SNR. The sky-average SNR is proportional 
to G and it is the same for each detector in the network 
(SNR w 2.3G). 

B. Simulation results 

In the simulation we tested the following detection 
methods: the standard likelihood method (L max ), the 
hard constraint method (ihard)j an d the soft constraint 
nethod (i so ft)- The detection performance of the meth- 
ods is compared by using the receiver operating charac- 
teristic (ROC), which shows the detection probability as 
a function of the false alarm probability. Examples of 
the ROC curves, corresponding to G = 3 and G = 4, are 
shown in Figure 

The accuracy of the source localization depends on the 
strength of the GW signal and on the network configu- 
ration. With only one detector in the network, the likeli- 
hood statistics is constant across the sky (it has no 6 or (f) 
dependence) and therefore the source localization is not 
possible. However, already with two spatially separated 
detectors the network becomes sensitive to the source lo- 
cation (see Figure |SJ| . In the case of two closely aligned 
detectors Hl-Ll, the area with the large values of the 
likelihood is rather a ring then a point, showing an am- 
biguity in the determination of the source location. But 
even in this case the method gives directional informa- 
tion about the source and allows exclusion of the most of 
the sky area as inconsistent with the detected GW signal. 
For two misaligned detectors Hl-Gl, the source localiza- 
tion is more accurate due to different angular sensitivities 
of the detectors. Even more accurate estimation of the 
source coordinates can be obtained with three and more 
detectors in the network (see Figure |3J). The greater the 
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FIG. 4: Receiver operation characteristics for the network of H1-L1-G1 detectors: SNR = 6.9 (left) and SNR = 9.2 (right). 
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FIG. 5: Sky maps of network statistics for 2 detector networks 
Hl-Ll (top) and Hl-Gl (bottom). The source is located at 
9 = 50° and <b = 280° . 



which satisfy the condition a < a c . The ratio of N a 
to the total number of injections defines the efficiency of 
the source localization and depends on the signal-to-noise 
ratio SNR. 

FigureHJshows the efficiencies of the source localization 
for the L1-H1-G1 network corresponding to the different 
detection methods. In this example, the values of the 
acceptable localization error are chosen to be a — 8° 
and a = 16°. Note that the constraint likelihood meth- 
ods perform considerably better than the standard like- 
lihood method. Let us consider, for example, the source 
localization for SNR of 10 (20) . The hard constraint 
method recovers approximately 48% (66%) of all simu- 
lated sources within the 8-degree angle from their true 
location. In comparison, the standard likelihood method 
yields only 12% (35%) efficiency for the same angle. 
Within the 16-degree angle, the hard constraint method 
recovers 66% (86%) of all the simulated sources, whereas 
the standard method yields only 22% (50%) efficiency. 
Similar comparisons hold for the soft constraint method. 
We find that for both constraint likelihood methods, the 
events with poorly reconstructed coordinates come from 
the areas in the sky with small values of the network 
sensitivity. 



number of the detectors in the network, the better the 
source localization. 

The error in the source localization is given by the 
angle a between the true direction to the source and the 
reconstructed direction to the source. Equivalcntly, a 
can be defined as the length of an arc connecting these 
two locations on a sphere with unit radius. To describe 
the efficiency of the source localization we introduce the 
following figure of merit. First, we chose a cone with the 
opening angle a c which constitutes an acceptable error. 
Then we calculate the number of detected sources (N a ) 



VIII. CONCLUSION 

We have presented a novel approach to the detection 
and reconstruction of gravitational waves with an arbi- 
trary network of interferometric detectors. Starting with 
the network likelihood ratio functional for unknown grav- 
itational wave burst signals, we identify and solve the 
two detector paradox. The essence of the paradox is 
that in the case of two arbitrary misaligned detectors 
the maximum likelihood ratio statistics depends only on 
the power in the detector data streams. It does not agree 
with the statistic of two co-aligned detectors, which de- 




signal-to-noise ratio signal-to-noise ratio 

FIG. 6: Efficiency of the source localization with a network of three detectors H1-L1-G1: a c — 8 (left), a c = 16 (right). 



pends also on the cross-correlation between the detec- 
tors. We show that the problem is associated with the 
different sensitivity of the detector network to two po- 
larization components of the GW signal and present not 
only in the case of two detectors, but for any arbitrary 
network. To characterize the difference in the sensitivity 
to the GW components, we introduce the network align- 
ment factor. For locations on the sky where the value of 
the alignment factor is small, the network is sensitive to 
only one GW component and the variation of the like- 
lihood functional results in the un-physical solutions for 
the second GW component. To exclude the un-physical 
solutions we propose to use constraints, which limit the 
parameter space of the GW waveforms and result in a 
new class of the maximum likelihood ratio statistics. For 
the networks of two and more detectors, the constraint 
likelihood methods allow reconstruction of the two GW 
polarization components and the location of the source 
on the sky. 

In the paper we introduce two examples of the con- 
straint statistics, which performance is compared with 
the standard likelihood statistics. The performance of 
the method was estimated with the numerical simula- 
tion. We restricted our simulation to the case of the 
white Gaussian noise. For simplisity we assumed that all 



detectors have identical sensitivities though the method 
presented in this paper does not have these restrictions. 
Our simulation results indicate that the constraint like- 
lihood method enhance the detection of the GW signals 
and performs significantly better than the standard like- 
lihood method in the reconstruction of the source coor- 
dinates. We believe that since all the methods we have 
considered are compared on exactly the same footing, our 
results regarding relative performance will not change for 
the general case. However, as a work in progress, we 
plan to expand our simulations to more realistic detector 
noise. 
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